load c:\manoj\projects\ace\pred_julia pred_w available_ace_fday missing_data_ace ace_matrix ;
load c:\manoj\projects\ace\Julia_W_new w_climate Julia_W  julia_model_drift fday;    
load c:\manoj\geomag\indices\aplist.mat;
%pred_w = pred_w.*24.366*1e-3; % electric fields in mV/m
w_climate = w_climate.*24.366*1e-3;
ace_matrix = ace_matrix./15;

number_of_days = 0;
%size of fday is 1002, number of ace overlap is 993 - this is because ace
%does not have data on the last 9 days of julia fdays 
h=figure(1);
set(h,'Position',[1 31 1920 1096]);
for i = 1:993,
    
        if length(Julia_W(i).k) >= 1,
            L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
            mean_ap = mean(ap(L));
        else,
            mean_ap=-1;
        end;

    if missing_data_ace(i) == 0 & available_ace_fday(i) == Julia_W(i).fday...
            & mean_ap >= 0,
      number_of_days = number_of_days+1;   
         obs_minus_climate = w_climate(i,Julia_W(i).k);
         obs_minus_climate = obs_minus_climate-nanmean(obs_minus_climate);
         obs_minus_climate_minus_ace = pred_w(i,Julia_W(i).k);
         obs_minus_climate_minus_ace = obs_minus_climate_minus_ace-nanmean(obs_minus_climate_minus_ace);
         climate_minus_ace = obs_minus_climate-obs_minus_climate_minus_ace;
         
         ace = ace_matrix(i,:);
         rms_error_1(i) = sqrt(mean(obs_minus_climate.^2));
         rms_error_2(i) = sqrt(mean(climate_minus_ace.^2));

win_dow =  rem(number_of_days,6);
    if win_dow==0,
        win_dow=6;
    end;         

subplot(3,2,win_dow);        
time_ax    = fday(i,Julia_W(i).k);     
plot(time_ax+datenum(2000,1,1),obs_minus_climate,'r','linewidth',1);
set(gca,'FontSize',10);
hold on;
%plot(time_ax+datenum(2000,1,1),obs_minus_climate_minus_ace,'b.');
plot(fday(i,Julia_W(i).k(1):Julia_W(i).k(end))+datenum(2000,1,1),...
    pred_w(i,Julia_W(i).k(1):Julia_W(i).k(end)),'b','linewidth',1);
plot(fday(i,:)+datenum(2000,1,1),ace,'k');
date_string = datestr(time_ax(1)+datenum(2000,1,1),29);
grid;
datetick;
axis([datenum(2000,1,1)+Julia_W(i).fday+(8+5)/24 datenum(2000,1,1)+Julia_W(i).fday+(18+5)/24 -inf inf]);
L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
mean_ap = mean(ap(L));
title(sprintf('Date - %s, Ap = %4.1f',date_string,mean_ap));
hold off;
if mean_ap > 40,
    pause;
end;

%pause;
if win_dow == 6| i==993,
% h=legend('Observed JULIA EF (mV/m)','Predicted JULIA EF (mV/m)','ACE Ey (mV/m)/15','Location','Best');
% set(h,'FontSize',8);
% saveas(gcf,['c:\manoj\projects\ace\plots\' date_string],'pdf');
% display(['c:\manoj\projects\ace\plots\' date_string]);
clf;
end;
end;
end;

a=nansum(rms_error_1);
b=nansum(rms_error_2);
fprintf('Number of days = %d, Percentage improvement = %10.5f\n', number_of_days,100*(a-b)/a);